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Abstract 

A  new  Discontinuous  Galerkin  Formulation  is  introduced  for  the  elliptic  reaction- 
diffusion  problem  that  incorporates  local  second  order  distributional  derivatives. 
The  corresponding  bilinear  form  satisfies  both  coercivity  and  continuity  properties 
on  the  broken  Hilbert  space  of  H2  functions.  For  piecewise  polynomial  approxima¬ 
tions  of  degree  p  >  2,  optimal  uniform  h  and  p  convergence  rates  are  obtained  in 
the  broken  H 1  and  H2  norms.  Convergence  in  L2  is  optimal  for  p  >  3,  if  the  compu¬ 
tational  mesh  is  strictly  rectangular.  If  the  mesh  consists  of  skewed  elements,  then 
optimal  convergence  is  only  obtained  if  the  corner  angles  satisfy  a  given  regularity 
condition.  For  p  =  2,  only  suboptimal  h  convergence  rates  in  L2  are  obtained  and 
for  linear  polynomial  approximations  the  method  does  not  converge. 
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1  Introduction 

Several  variations  of  Discontinuous  Galerkin  Methods  (DGM’s)  for  second  or¬ 
der  elliptic  boundary  value  problems  have  been  proposed  during  recent  years, 
which  exhibit  special  convergence,  conservation  and  local  approximation  prop¬ 
erties  attractive  for  parallel  adaptive  /ip-approximations.  A  comprehensive  ac¬ 
count  of  several  types  of  DGM’s  can  be  found  in  the  volume  edited  by  Cock- 
burn,  Karniadakis  and  Shu  [1],  in  the  paper  of  Arnold,  Brezzi,  Cockburn,  and 
Marini  [2],  and  in  the  report  by  Prudhomme  et  al.  [3]. 

In  1997,  Oden,  Babuska,  and  Baumann  [4]  introduced  a  discontinuous  Galerkin 
formulation  similar  to  the  GEM  formulation  by  Delves  et  al.  [5,6],  but  sign 
differences  in  certain  terms  resulted  in  a  positive  definite  bilinear  form.  For 
p  >  2  (where  p  denotes  the  minimum  order  of  the  polynomial  approximations), 
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the  DGM  by  Oden,  Babuska,  and  Baumann  appears  to  be  unconditionally 
stable,  whereas  the  GEM  formulation  requires  the  inclusion  of  penalty  terms 
to  stabilize  the  formulation.  Complete  details  are  given  in  [7]. 


Riviere  et  al.  [8-10]  proposed  an  extension  of  the  DGM  of  Oden,  Babuska,  and 
Baumann  by  including  a  penalty  term  on  the  jumps  of  the  solution  across  the 
element  interfaces.  The  formulation,  due  to  the  addition  of  a  penalty  term, 
does  not  satisfy  a  local  conservation  property  any  longer,  but  it  becomes  stable 
for  p  =  1.  Moreover,  optimal  h  and  suboptimal  p  convergence  rates  are  proved 
in  the  broken  H 1  space. 


In  this  paper,  a  new  DGM  formulation  is  presented  for  a  model  class  of  linear 
elliptic  boundary  value  problems,  that  incorporates  the  second  order  distri¬ 
butional  derivatives,  and  consequently  is  defined  within  the  space  of  local  H 2 
functions.  The  formulation  does  not  require  any  penalization  or  stabilization, 
exhibits  a  local  conservation  property,  and  guarantees  optimal  uniform  h  and 
p  convergence  in  the  broken  H 2  and  H 1  spaces  when  the  order  of  polyno¬ 
mial  approximation  is  at  least  of  degree  >  2.  This  is  an  improvement  over 
other  DGM  formulations  for  which  optimal  h  convergence  has  been  proved 
but  optimal  p  convergence  rates  have  not  been  established. 


This  formulation  exhibits  eccentric  convergence  behavior  in  L2 .  When  the 
computational  meshes  employ  solely  rectangular  elements,  it  is  proved  (and 
numerically  confirmed)  that  optimal  h  and  p  convergence  rates  hold  for  poly¬ 
nomial  approximations  of  degree  >  3.  For  skewed  meshes,  we  can  prove  that 
under  certain  conditions  the  error  can  converge  optimally.  However,  in  practi¬ 
cal  computational  applications,  we  expect  such  cases  to  be  rare  and  one  should 
generally  expect  suboptimal  convergence  rates  in  L2(D). 


This  is  not  the  first  DGM  that  uses  second  order  derivatives.  In  the  works  by 
Engel  et  al.  [11,12],  these  derivatives  are  also  incorporated  by  adding  Galerkin 
Least  Squares  (GLS)  stabilization  terms  to  existing  DGM’s.  The  DGM  in¬ 
troduced  in  this  paper  starts  with  GLS  terms  on  each  element  that  are  the 
variational  equivalent  of  a  model  reaction-diffusion  problem.  The  final  for¬ 
mulation  is  then  derived  by  applying  classical  Green’s  identities  (see  proof 
of  Lemma  2)  to  the  GLS  terms  and  enforcing  continuity  across  the  element 
interfaces  and  boundary  conditions  on  the  outer  boundary. 


The  outline  of  this  paper  is  as  follows:  in  Section  2,  the  model  problem,  no¬ 
tations,  the  proper  function  space  setting,  and  the  new  DGM  are  introduced. 
Subsequently,  the  convergence  properties  of  the  DGM  are  proved  in  Section  3. 
In  Section  4,  the  theoretical  results  are  confirmed  on  one-  and  two-dimensional 
problems.  Section  5  concludes  with  a  brief  summary  of  results. 
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2  The  New  DG  Formulation 

2.1  Model  Problem  and  Notations 

Let  C  M2  be  an  open  bounded  domain  with  Lipschitz  boundary  dfl,  and 


—  rD u rN,  rD nr^v  —  0 , 


where  T o  denotes  the  part  of  the  boundary  with  prescribed  Dirichlet  boundary 
conditions  and  T n  the  part  subjected  to  flux,  or  Neumann,  conditions.  Let  Vh 
be  a  regular  partition  of  hi  into  N  open  elements  {A",}  with  diameters  {hi}, 
such  that  (see  Figure  1): 


(1) 


The  maximum  diameter  in  the  partition  is  denoted  as  h  and  the  unit  normal 
vector  n*  on  the  edge  dKj  of  Ki  G  Vh  is  directed  outward  with  respect  to 
the  element  Kt  (see  Figure  1).  Given  an  element  Kj  G  Vh,  the  part  of  the 
boundary  of  dKj  that  is  shared  with  a  neighboring  element  Kj  is  denoted 
dKij,  i.e. 


(2) 


dK^  =  dKi  n  dKj. 


Note  that  dK y  is  an  open  subset  of  dKj .  We  now  consider  the  following 
reaction-diffusion-type  model  problem: 


Find  u,  such  that  : 


—A u  +  u  —  /,  in  Q 


(3) 


u  =  Wo,  c»n  r £> 
V«  •  n  =  g,  on  T at 


Fig.  1.  Geometrical  definitions 
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where  /,  the  source  term,  is  a  real- valued  function  in  A2(h2),  and  the  Dirichlet 
and  Neumann  boundary  data,  u0  and  g,  are  respectively  in  Hi^2{dKi  D  Tp) 
and  in  Hl^2(dKi  D  Tjv),  for  all  K%  G  Vh  whose  edges  intersect  with  dVt. 


2.2  The  Broken  Banach  Space  of  Test  Functions 

Given  the  partition  Vh,  the  following  broken  space  is  defined: 


H2(Vh )  =  ju  G  A2(fl)  :  v\Ki  G  H2(Ki ),  VAT*  G  A,  j, 


(4) 


We  introduce  two  norms  on  this  space, 


HhHVu)  =  E{llV2'‘’lli.«)  +  2  \\Vv\\lHKl)  +  ||Hll.(K)J, 

iib'ii2  —  y  ( + 2  iivvii^2(K.)  +  imi2^)}. 


i=  1 


(5) 


where  the  first  uses  the  local  Sobolev  norms  in  H2(Ki )  and  the  latter  the 
Laplacian  norms.  Within  this  setting,  we  introduce  the  local  zeroth  and  first 
order  trace  operators  [13,14]  for  functions  on  K,  G  Vh' 


To  : 

H\K[)  - 

-»•  Hl/2{dKi 

7i  : 

H2{ Ki)  - 

-»•  H^2(dKi 

(6) 


where  (vf)  represents  the  trace  of  the  normal  derivative  dvi/dn  on  dK% .  We 
define  the  norm  on  Hl^2{dKi)  as  follows  [13]: 


IMI  HV*idKi)  =  IMUw, 

7q  («)=¥> 


Remark  1  (Trace  Inequalities)  Let  Vh  consist  of  elements  {Ki}  with  Lips- 
chitz  boundaries.  Then,  the  trace  operators  (6)  are  continuous  and  surjective 
(see  [13]), 


WYo^Wh1/2 (dKi)  <  IKIIrP(iy), 
Il7l(v*)ll HWidKi)  <  C  \\vi\\H2(Ki),  C  >  0. 


(7) 
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2.3  The  Weak  Formulation 


The  discontinuous  variational  formulation  is  stated  as  follows: 


Find  w  G  H2(Vh)  : 

B(w,  v)  =  L(v),  Vu  G  H2(Vh ), 


(8) 


where  the  bilinear  form  B(-,  •)  and  linear  form  L(-)  are  defined  as: 

B  :  H2(Vh)  x  H2{Vh)  — >  R,  L  :  if2(n)  — >  R, 


N 


i= 1 


IK;  L 


B{w ,  v)  =  E  /  AwiAvi  +  2  V«)j  •  Vt)j  +  WjUj 


dx 


E 


dKijCdKi  JdKH 


TiK'boOi)  -7iE)7oK) 


ds 


/  .  7i(l,i)  To(wi)  ds  —  /  _  7i(«k)7obi)  ds}, 

J  dKiC\Y]3  J 


(9) 


f(-Avi  +  Uj)  dx  + 


/a/v(nrD 


7i(^i)«o  ds 


+ 


/afTinrjv 


P7o(vi)  ds 
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where  denotes  U|xr  Analogous  to  the  DG  formulation  by  Oden,  Babuska, 
and  Baumann  [4],  this  formulation  satisfies  local  balance  of  conservation  laws. 
By  taking  a  function  v  such  that  v  —  1  on  an  element  Kt  G  Vh  and  v  =  0 
elsewhere,  (8)  and  (9)  give  us: 


idi<inrD 


ds 


'Kt 


/dx-  E 


dKijCdKi  JdKH 


7i  (wj)  ds 
g  ds. 


idKtn  rN 


Let  w  be  a  solution  of  (8).  Then  the  left  hand  side  of  the  above  expression 
represents  the  reaction  ‘forces’  of  the  element  to  external  body  and  boundary 
loads  (right-hand  side).  Moreover,  we  can  prove  that  w  satisfies  the  PDE  (3) 
and  belongs  to  the  following  space: 


if  (A,  ft) 


G  L2{VL)  :  Av  G  L2{Vl) 


Lemma  1  Let  f  G  L2(ft),  g  G  Hl/2(dKiC\T N)  andu0  G  H3^2(dKir\TD) ,  and 
suppose  w  G  H2(Vh)  is  a  solution  to  the  discontinuous  Variational  Boundary 
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Value  Problem  (VBVP)  (8).  Then,  w  satisfies  the  PDE  (3)  in  the  distribu¬ 
tional  sense  and  belongs  to  H( A,fl)  fl  H2(Vh). 

Proof:  By  taking  an  arbitrary  K,  G  Vh  and  choosing  a  smooth  test  function  v 
that  vanishes  outside  Kt,  then  (8)  and  (9)  yield: 


iKi  L 


AwiAi^i  +  2  Vwi  ■  Vpt  +  Wi(pi  dx  =  /  f(-Api  +  pf)  dx. 


IKi 


Application  of  Green’s  first  identity  to  the  integral  of  2  'Vwl  ■  'Vifi,  leads  to: 


{-Awi  +  Wi) 


{-A ^  +  pi)  dx 


f(~Api  +  dx,  Wpi  G  V(Ki). 


Thus,  in  the  distributional  sense,  w  satisfies  the  following  problem  on  any 
Ki  G  Vh: 


-A  Wi  +  Wi  =  f. 


(10) 


We  return  to  (8)  and  (9)  and  again  consider  an  arbitrary  element  K,  G  Vh 
and  the  interface  dKtJ  of  this  element  with  one  of  its  neighbors  Kj.  We  choose 
test  functions  p  that  vanish  outside  of  Ki  and  everywhere  on  dKl}  except  on 
the  element  interface  dKVJ .  Substituting  such  test  functions  and  then  applying 
Green’s  first  identity,  yields: 


/ Ki 


(-A w  +  w)  (-A pi  +  pt)  dx  +  /  ^  Wiiwi)  +7i(wt)l  7o(<A)ds 


dKi 


+  /  7i  (<Pi)  7o(wi)  -  7o(wj)  ds  =  f(—Api  +  p^  dx. 


I  OK, 


Ki 


By  recalling  (10),  this  expression  gives  the  weak  continuity  of  w  and  its  normal 
flux  dw/dn  across  the  element  interface  dK,p. 


7o  (wi)  =  7o  (wj ) )  7 lOi)  =  -7i  (wj). 


Obviously,  by  repeating  this  procedure  we  establish  the  weak  continuity  of 
w  and  dw/dn  across  any  element  interface  dKij  in  the  partition  Vh,  which 
implies  that  w  is  in  H(A,Q).  To  prove  satisfaction  of  the  Neumann  boundary 
condition,  we  take  test  functions  p  that  vanish  outside  of  K \  and  everywhere 
on  dKi  except  on  that  part  of  dK. \  that  coincides  with  the  boundary  T jy.  By 
using  such  test  functions,  we  now  get: 


(-Aw  +  w)  (—Api  +  p^  dx  + 


ldKinrN 


7i(wi)  7o(Pi)  ds 


=  /  / (~Api  +  Pi)  dx,  +  /  97o(Pi)  ds. 

J  Ki  JdKi  nrN 
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Again,  recalling  (10)  reveals  that  the  Neumann  boundary  condition  on  dK,  n 
T jv  is  satished  weakly: 


/  7l(^i)7o(^)ds  =  /  #7o(<A)ds. 

/  dKif)T  jy  J  dKiDT  jy 


Analogously,  we  can  prove  satisfaction  of  the  Dirichlet  condition. 


In  the  next  lemma,  we  prove  the  converse  of  Lemma  1  and,  therefore,  establish 
equivalence  between  the  weak  and  strong  formulation  of  the  model  problem. 

Lemma  2  Let  u  G  H(A,Q)  fl  H2(Vh )  be  the  solution  to  problem  (3),  then  u 
is  a  solution  to  the  VBVP  (8). 

Proof:  By  taking  an  arbitrary  test  function  G  H2(Ki),  multiplying  (3)  by 
(—A Vi  +  v^,  and  integrating  over  K gives  us  the  Galerkin  Least  Squares 
(GLS)  representation  of  (3)  on  K,: 


/  ,  (-A m  +  Ui)  (-A m  +  Vi)  dx  =  f  (-A n  +  vf)  dx. 

Summing  the  contributions  for  all  elements  in  Vh,  yields: 

N  f  N  f 

X!  /  ,  (-^Ui  +  Ui)  (-A Vi  +  Vi)  dx  =  /  /  (-An*  +  Vi)  dx. 

i=\  i=  1 

By  applying  Green’s  first  identity  to  the  integrals  of  UiAvi  and  A UiVi,  we  get: 
/  [AuiAvi  +  2Vui  -Vvi  +  u^i }  dx 


i=i 


f  [l\(Ui)lo(Vi)  +  Yl(Vi)Y0(Ui) 
JoKi 


N 


ds  f  =  /  ,  f  (-A Vi  +  Vi)  dx. 

J  i=i  -'A 

HD 

Concentrating  on  the  edge  integrals  and  applying  the  boundary  conditions, 


7o (ui)  =  u0,  on  dlu  nrD, 
7i  (ui)  =  g,  on  dlu  n  Tn, 


leads  to: 


IdKi  L 


7i(a)7o(a)  +  7i(a)7o(a 


ds=  /  7i(«i)7o(A)ds 

JdKi\TN 


(12) 


+  /  71(a)  7o(«i)ds  +  /  SAo (a)  ds+  /  71(a)  w0  ds. 

JdKi\rD  JdKiCTN  JdKinTo 
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Since 


>dKi\ rN 


Yi(ui)  Yo(vi)  ds 


I]  7l(«i)7o(y*)ds 

dKijCdKi  JdKii 


IdKiDFa 


7i(«i)  7o(u»)  ds, 


ldKi\Ti 


ll{vi)  lo(ui)  ds 


5]  L  7i(u<)7o(w<)ds 

dKijCdKi  JdKa 


JdKirTN 


7i(^)7o(M*)ds, 


we  can  rewrite  the  right-hand-side  of  (12),  which  gives: 


IdKi  L 


7lKbo(^)  +7Kut)75(«i) 


ds  = 


E 


dKijCdKi  JdKH 


7i  («*)  7o(yi)  +  7i(«i)  7o(«i)  ds  +  /  7i(w<)  7o(w;)  ds 


+  /  7i(v*)  7o('u*)  ds  +  /  07o(uO  ds  +  /  7^)  w0  ds. 

JSA'inr^  iaAinrjv  2aA'inrD 

(13) 

The  solution  u  to  (3)  is  in  H(A,  O)  and,  therefore,  its  trace  and  normal  deriva¬ 
tives  across  the  element  interfaces  dK^  are  continuous,  i.e. 


7o  (M* )  =  7o  ( ) ,  7i  («i )  =  -7i  («j )  r  on  every  dKi:j , 


where  Uj  denotes  the  restriction  of  u  to  neighboring  Kj,  and  7 l(uj)  and  7) (uj) 
are  the  traces  of  Uj  on  dKj.  Implementing  these  continuity  conditions  in  (13), 
yields: 


IdKi  L 


fl(Ui)Y0(Vi)  +  Yl(Vi)f0(Ui) 


ds  = 


E 


dKijCdKi  *  dKij 


1l(uj)  7o(Vi)  ~  7l(vi)  7o(uj) 


ds+  /  7i(ui)  7o(vi)  ds 
JdKinrD 


/  lfi(vi)  lfo(ui)  ds  +  /  gYoivi )  ds  +  /  7i(ui)uods. 

/9Ainrw  Jmynrjv  JdKirTD 


By  substituting  this  result  back  into  (11)  and  by  noting  that 


N 

E  E 


N 


%<=  1  dKijCdKi  dKii 


li{vi)  7o(Mi)  ds  =  E  E  L  7i  E)  7o('u*)  ds, 


i=l  dKijCdKi  dKii 


we  finally  get: 


N 


U  L  [AuiAvi  +  2 'Vui  ■  'Vvi  +  UiVi]  dx 


i= 1 


IK 


E 


dKijCdKi  "dKii 


7l  (uj)  7o(g)  ~  7l  (vj)  Ho(ui 


ds 


f  71(g)  7o(g)  ds  -  /  71(g)  7o(g)  ds  1 

J  dKiHT  £)  J  dKiDT  j\j-  J 

=  E  {  /  /  (~AG  +  g)  dx  +  f  g  %(vi)  ds  +  /  7i  (g)  uo  ds)  . 
1  JdKinrN  JdKinrD  J 


i=  1 


IKi 


which  establishes  the  assertion.  ■ 

The  proof  of  Lemma  2  uses  a  GLS  representation  of  (3)  on  every  element  K, . 
Engel  et  al.  [11,12]  have  used  such  GLS  terms  in  their  DGM  formulations  as 
a  stabilization  to  classical  DGM’s  (e.g.  [4,9,15,5]).  Here,  however,  they  serve 
as  a  starting  point  in  the  derivation  of  a  DGM  formulation,  whose  final  form 
is  obtained  by  applying  Green’s  identities  to  enforce  continuity  and  boundary 
conditions. 


2-4  Continuity  and  Coercivity  Properties 

The  functionals  and  L(-)  satisfy  continuity  and  coercivity  properties 

on  the  space  H2(Vh)-  We  start  with  an  important  coercivity  property  of  the 
bilinear  form  !?(•,•)  on  H2(Vh)  x  H'2(Vh)- 

Lemma  3  Let  -£>(•,•)  be  the  bilinear  form  defined  in  (9).  Then,  is  co¬ 

ercive  with  respect  to  the  broken  Laplacian  norm  |||  •  ||| , 

B(v,v)  >  ^l\v\\2,  VveH2{Vh). 


Proof:  Taking  w  —  v  in  (9),  yields: 


N 


B(v,v)  —  El  IIAd||!2(^)  +  2  ||Vnj|||2(X.)  +  IKIlE/f.) 


2—1 


ldKiD(TNurD) 


71(g)  75(g)  ds 
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Since  Yo(vi)  e  H^2(dKi)  and  Yivi  G  H  1^2{dKi)  for  all  K%  G  Vh,  the  Cauchy- 
Schwarz  Inequality  can  be  applied  to  the  boundary  integrals,  leading  to: 

B(v,v)  >  51]  \\Avi\\2L2(K.)  +  2  ||  Vui||22(A'.)  +  \\vi\\'2L2(K.) 
i= 1 1 

—  \\ll{Vi)\\H-1/2(dKi)  \\Yo(Vi)\\H1/2(dKi)i  }■ 


Applying  Young’s  inequality,  gives: 

N  ( 


B(v,v)  >  E]  \\AVi\\h(Ki)  +  2  llVuiUi^x-.)  +  \\vi\\2L2{Kz 


i— 1 


'9  1171(^)11^-1/2(9^)  9  ll7o(v*)ll//i/2(9iCi) 


Recalling  the  trace  inequality  (7)1,  yields: 


N 


B(v,v )  >  EjllA^H!2^)  +  2  IEv*lli2(A'i)  +  2  11^*11  h2(A'i) 


9  ll7l('yi)llE1/2(9A'i) 


(14) 


We  call  upon  Theorem  2.5  in  [13]  and  state  the  following  trace  inequality  for 
v  G  H2(Kty. 

1171(^)11^-1/2(9^)  <  II Auilll2^.)  +  ||Vu,||!w.  (15) 

Substituting  this  inequality  into  (14)  establishes  the  assertion.  ■ 


We  cannot  prove  coercivity  of  the  bilinear  form  with  respect  to  the  broken 
Sobolev  norm  IHIr/2^)-  With  the  issue  of  continuity,  we  face  a  converse  sit¬ 
uation:  continuity  of  B( •,  •)  on  H2(Vh)  x  H2(Vh)  can  be  proved  in  the  norm 
ll'll//2(ph)  but  not  in  the  norm  |||-|||. 

Lemma  4  The  bilinear  form  B(-,-)  is  continuous  on  H2{Vh)  x  H2{Vh),  he. 
there  exists  a  constant  M  >  0  such  that: 


\B(w,v)\  <  M\\w\\H2{Vh)  \\v\\H2{Vh)  Ww,v  G  H2{Vh). 


Proof:  Since  w,  v  G  H2(Vh),  the  zeroth  and  first  order  traces  of  these  functions 
are  respectively  in  Hi^2{dKf)  and  which  are  both  subspaces  of 
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L2(dKi).  Thus,  we  can  apply  the  Cauchy-Schwarz  inequality  to  £>(•,  •)  in  the 
following  manner: 


B(w,v)  <  53|llA^ilU2(Ki}  \\AVi\\L2iKi)  +  2  \\Vwi\\L2{Ki)  ||V^||L2(^.) 

+lh»ll L2(Ki)  II^IIl2^)  +  E  {IIti^OIIl2^,)  I|7oK)IU2(9^) 

dKijCdKi 

+  W'yi(vj)\\L2{dKij)  ||7o(wi)IU  2(dKij)}  +  W'yl(vi)\\L^{dKinrN)  1 1 7o  ( )  I U2  (^nrjv ) 


+ 1 1 7i  ( )  1 1  l2  nrD )  1 1 7o  K )  I U2  nrD ) 

We  can  bound  this  inequality  as  follows: 


B(w,v)  <  C 


\ 


N 

hill2  +  E  {hi(wi)\\h (dKi)  +  H(wi)\\h{aKi)} 


i=  1 


X 


N 

\  IMP  +  E  {  Il7l  (^i)  llz,2(6»JftTi)  A  ll7o('ui)llz,2(,9xi)}’  C  >  0, 

\  i=  1 


where  |||'j|j  is  defined  as  in  (5).  Since  H1B(dKi)  is  embedded  in  L2(dKj),  we 
can  assert  that: 


B(w,v)  <  C 


\ 


N 

hill2  +  E  { IN  hO  II  #  1/2^)  + 1 1 7^  K)  11^1/2^)} 


i=  1 


N 


X, 


Nil2  +  E  {ll7ih)ll^i/2(^.)  +  H(vi)\\2m/2{dKi)},  C>0- 

\  i=  1 

Application  of  the  trace  inequalities  (7)  completes  the  proof.  ■ 

Proposition  1  The  linear  form  L(-)  is  continuous  on  H2(Vh): 

30  0:  L(v)  <  C  \\v\\H2(Ph),  Wv  G  H2(Vh),  (16) 

where  C  =  C(f,  u0,g). 


Remark  2  (Well  Posedness)  Although  the  bilinear  form  B(-,  •)  satisfies  coer- 
civity  and  continuity  properties,  we  cannot  invoke  the  Generalized  Lax-Milgram 
Theorem  to  prove  existence  of  unique  solutions  of  (8),  for  two-dimensional 
problems.  The  coercivity  property  in  Lemma  3  is  satisfied  in  terms  of  the 
norm  |||-|||.  For  two-dimensional  problems,  the  space  H2{Vh )  is  not  complete 
with  respect  to  this  norm  and  completeness  is  an  essential  condition  in  the 
Lax  Milgram  Theorem. 
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However,  for  classes  of  problems  of  (3)  for  which  we  can  prove  there  exists  a 
solution  u  G  H2(Vh)  f~l  H(A,Q)  (e.g.  via  the  conventional  continuous  varia¬ 
tional  formulation) ,  we  know  from  Lemma  2  that  u  is  a  solution  to  the  DGM 
formulation  (8).  Uniqueness  is  then  guaranteed,  as  the  bilinear  form  is  positive 
definite  on  H2(Vh )■ 


3  Convergence 

3.1  The  Discrete  Problem 

Let  {F^}  be  a  family  of  invertible  maps  defined  on  the  partition  Vh  such  that 
every  element  Kt  G  Vh  is  the  image  of  F k,  acting  on  a  master  element  K,  as 
shown  in  Figure  2. 

Fa-,  :  K  — >  Ku  x  =  F*i(x).  (17) 

Unless  stated  otherwise,  the  sets  of  mappings  are  assumed  to  be  affine.  We 
introduce  a  finite  dimensional  space  of  real- valued  piecewise  polynomial  func¬ 
tions, 

V*"  =  {v  e  L2(Sl)  :  =  v  o  F<y  ,  v  e  PP‘(K),  V/t,  e  Ph\  C  H2 (Vh), 

(18) 

where  PPl  (K)  denotes  the  space  of  polynomials  on  K  of  degree  <  jy ,  in  which 
Pi  can  have  different  values  on  different  elements.  Let  u  G  H2(fPh)  be  the 
solution  of  (8).  Then,  we  seek  a  discrete  approximation  Uh  G  Vhp  by  solving 


Fig.  2.  Mapping  from  the  master  elements  to  the  physical  space 
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the  following  (discrete)  variational  problem: 


Find  uh  G  Vhp  : 

B(uh,vb)  =  L(vh),  Vm  e  v'lp. 


(19) 


Lemma  5  The  bilinear  form  £?(-,•)  is  coercive  on  Vhp  x  Vhp  with  respect  to 
the  norm  ||-||jj2(ph)  (see  (5)),  i.e. 


30  0:  \B(vh,vh)\  >  C  \\vh\\2H2{Ph),  Vvh  G  Vhp. 


Proof:  Since  Vhp  C  H2(Vh),  we  know  from  Lemma  3  that  the  bilinear  form 
B(-,  •)  is  coercive  on  Vhp  x  Vhp  with  respect  to  the  norm  |||-|||.  For  finite  dimen¬ 
sional  spaces,  this  norm  is  equivalent  to  |M|H2(pfe)-  Thus, 

3C  >  0  :  |K||»2(J>„)  <  c||t>i|||,  Vui  e  vhr, 

which  establishes  the  assertion.  ■ 


Since  the  bilinear  form  B(-,-)  is  continuous,  coercive,  and  positive  definite 
on  Vhp  x  Vhp,  existence  of  unique  solutions  Uh  G  Vhp  to  (19)  is  established 
by  applying  the  Generalized  Lax  Milgram  Theorem.  If  u  is  then  the  solution 
to  (8),  it  easily  follows  that  the  approximation  error  =  u  —  is  governed 
by  the  following  variational  problem: 


Find  eh  G  H2(Vh)  such  that 

B(eh,v)  =  L(v)  -  B(uh,v ),  Vv  G  H2(Vh ) 

Tlhiy) 


(20) 


where  IZh  '■  H2(Vh)  — >  M  is  the  Residual  Functional,  which  satisfies  the 
Galerkin  orthogonality  property  on  the  space  Vhp,  i.e. 


B(eh,vh)  =  Hh(v)  =  0,  \tvh  G  Vhp. 


(21) 


3.2  A  Priori  Error  Estimates  in  Fd2{Vh) 

In  this  section,  we  derive  convergence  rates  of  the  approximation  error  eh  = 
u  —  Uh  in  terms  of  the  norm  1 1  - 1 1  ^ .  The  convergence  rates  in  lower  norms 
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are  derived  in  the  next  section.  We  start  by  defining  a  set  of  interpolants  {p lhp} 
for  every  Vh,  such  that: 

nip  :  Hri(Ki )  — »  PPi(Ki ),  K,  G  Vh,  i  =  1,  2, . . . ,  N, 

Kp(vh )  =  Vh,  V vh  G  PPi ( Ki ) , 

where  r,  >  2.  We  can  now  call  upon  an  interpolation  theorem  proved  in  [16]. 

Theorem  2  For  p>  G  Hn(Ki),  there  exists  C  >  0,  independent  of  (p,  p:  and 
Ti,  and  a  sequence  n lhp{(p)  G  PPi(Ki),  such  that: 


By  extending  the  local  interpolants  Fhp(.)  to  zero  outside  of  Afi ,  we  can  define 
a  global  interpolant  on  the  whole  partition  TV 


nv :  H\vh)  ~  Vh“,  nv(t.)  @  •£  Oik).  » e  (22) 


Lemma  6  (Interpolation  Lemma)  Let  u  G  H2(Vh).  Then,  there  exists  C  >  07 
independent  ofu,  {hi}  and  { pt }  such  that  the  interpolation  error  77  =  u  —  II 
can  be  bounded  as  follows: 


INI*™  < 


X!  IMIfrqjqp.  ri  P  2, 


where  r 


min  {rj},  p  =  max  {pi},  h  =  max  {/q},  and  /i  =  min(p  +  l,  r 
KieVh  h 


Proof:  The  proof  of  this  lemma  is  quickly  established  by  recalling  the  definition 
of  the  norm  1 1 •  1 1 7/2 (vh ) ,  as  given  in  (5),  and  substituting  the  inequalities  listed 
in  Theorem  2.  ■ 


Having  established  convergence  rates  for  the  interpolation  error,  we  can  now 
derive  optimal  convergence  rates  of  the  approximation  error  in  the  broken 
space  H'2(Vh). 
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Theorem  3  ( Convergence )  Let  u  G  H2(fPh)  be  the  unique  solution  to  the 
variational  problem  (8)  and  {uh  G  Vhp}  be  a  sequence  of  approximations  (19) 
of  u.  Then,  the  approximation  error  en  =  u  —  Uh  is  bounded  as  follows: 


^  ^ 

eh\\H2(Ph)  <  C  2  /  IMI/Fi(iq)J  Pi  —  1  >  ri  —  2, 

p  V  Ki^h 


where  r  =  min  in),  p  =  max  \pi\,  h  =  max  {hi],  and  u  =  min(p  +  l,r). 

K^vhXli,F  iu&vhXFlS’  K^vhy  ^  J 

Proof:  Let  n^p(-)  be  the  interpolant  operator  as  defined  in  (22).  Then  we 
introduce  two  functions  r)  and  £,  such  that  the  approximation  error  can  be 
written  as  =  rj  —  £,  where  77  =  u  —  HhPu  and  f  =  Uh  —  II hpu.  Note  that 
f  g  Vhp  and  that  the  interpolation  error  r]  is  in  H2(Vh)-  By  using  the  triangle 
inequality,  we  obtain: 


\\eh\\H2(vh)  <  IMI H2(rh)  +  Ill'll H2(vh)-  (23) 

Recalling  the  coercivity  property  of  Lemma  5  leads  to: 

1  m2H^vh)<cB^i). 

By  applying  the  orthogonality  property  (21)  and  the  continuity  of  the  bilinear 
form  B(-,  •)  (see  Lemma  4),  we  get: 

Ill'll H2(vh)  <  CB(r),£)  <  C  \\r)\\H2(vh)\\€\\H2(vh) 

Thus,  returning  to  (23),  we  can  conclude: 

\\eh\\H2(vh)  <  C\\v\\ H2(vh)- 

We  finish  the  proof  by  applying  Lemma  6.  ■ 

Remark  3  Take  p  =  1.  From  Theorem  3,  p  =  2  and  /i  —  2  =  0,  which  implies 
that  the  approximate  solutions  {uh}  do  not  converge  for  h  refinements. 


3.3  The  Aubin-Nitsche  Lift  -  A  Priori  Error  Estimates  in  Lower  Norms 

For  p  >  2,  convergence  to  the  solution  of  the  target  problem  (8)  is  guaranteed 
by  Theorem  3,  but  the  rates  are  suboptimal  in  terms  of  the  E[1( Vh)  and  L2{Et) 
norms.  We  employ  a  technique  introduced  by  Aubin  [17]  and  Nitsche  [18]  to 
prove  that,  under  certain  conditions,  the  approximation  error  also  converges 
optimally  in  these  lower  norms.  First,  we  introduce  broken  Hilbert  spaces  on 
the  polygon  partitions  {Vh}- 

H°{Vh )  =  {v  G  L2(Ll)  :  v\Ki  G  VAfi  G  Ph)  ,  0  <  a  <  1,  (24) 
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on  which  we  define  broken  Sobolev  norms, 


Ml  H°{Vh) 


\ 


£  I 

Ki&Vh 


12 
I  Hc 


{KiV 


e  E°{Th 


(25) 


We  follow  [17,18]  by  introducing  the  functionals  q(-)  in  the  dual  space  H~a(Th ), 
for  which  we  can  prove  that  there  exist  Riesz-type  representative  functions  wq 
in  the  following  subspace  of  H2(fPh)' 


HUVh) 

Lemma  7 

that: 


=  ju  G  H2(Vh)  :  7oM) 
For  every  q  G  H~a(Vh) 


=  0,  7i(Wi)  =  0,  VA7  G  C  i/2(n). 
t/iere  exists  a  unique  wq  G  H2]{q(Vh),  such 


B(v,w,)  -  q(v),  Vv  e  H2{Vh) 


(26) 


Proof:  If  we  recall  the  coercivity  property  of  Lemma  3,  we  get: 


\B(v,w)\  1  I  w|||2 

sup  ~ nil -  —  o  — M il - 

v&H^{Vh)/{  0}  \\v\\H2(Vh)  2  \\w\\H2{Ph) 


Vu>  e  H2m(Vh)/{ 0}. 


For  functions  that  belong  to  H^Kf),  the  norm  \\Awl\\[2(K.j  is  equal  to  the 
norm  || V2Wj||^2(X-i)  ( e.g .  see  [13]).  From  (5)  then  follows  that  the  norms 
INI H2(vh)  and  ll'l  are  identical  for  functions  that  belong  to  HffVh)-  Hence, 
the  above  expression  gives  us  the  Inf-Sup  condition  for  the  bilinear  form  B( •,  •) 
on  H2(Vh)  x  H20(Vh): 


sup 

veH2(ph)/{0} 


B(v,w)  | 
IMI  H2{Ph) 


>  2  IHIff2(Ph)> 


Vw  G  H20(Vh)/{ 0}. 


Considering  that  the  bilinear  form  is  also  positive  definite  and  continuous  (see 
Lemma  4),  we  can  call  upon  the  Generalized  Lax  Milgram  Theorem  to  assert 
that  there  exists  a  unique  solution  wq  G  H2m(Vh).  ■ 


By  duality,  the  norm  of  the  error  in  the  spaces  Ha(Vh),  0  <  cr  <  1,  is  closely 
related  to  the  functions  wq,  i.e. 


¥h\\H°{Vh)  =  sup 

qeH—(ph) 


I  q(eh)\ 

Ikll  H-^(Vh) 


sup 

q&H-°(Th) 


\B(eh,wq)  | 

\\<l\ \h-°(vh) 


Let  nhp(-)  G  Vhp  denote  the  global  interpolation  operator  as  defined  in  (22). 
Then,  by  applying  the  Galerkin  orthogonality  property  (21)  and  continuity  of 
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the  bilinear  form,  we  can  rewrite  this  expression  as: 


\eh\\m{vh)  <  C  \\eh\\H2(vh)  sup 

qdH-°{Vh) 


u>n 


Ikll  H-°(Vh) 


C>  0. 


By  applying  the  interpolation  Lemma  6  and  convergence  Theorem  3,  we  can 
further  bound  the  error, 


\\e-h\\H°{vh)  -  C 


v — 4 


.r+s— 4 


p 


/  IMIlrqiq) 

KidVh 


sup 

q£H-°(Vh) 


I  \\wq\\2Hsi(Ki) 

Ki&vh  _  (27) 

H-°(Vh)  ’ 


where  h ,  p,  r,  r*,  and  p  are  dehned  in  Theorem  6,  and  s  =  min  {sj}  and 

KidVh 

v  =  min(p  +  l,s).  So  convergence  of  the  approximation  error  in  the  lower 
Hl(Vh)  and  L2(Q)  norms  is  governed  by  the  regularity  of  the  solutions  wq.  To 
determine  this  regularity,  we  call  upon  a  regularity  theorem  that  is  based  on 
the  work  on  polygonal  domains  by  Grisvard  [19]. 

Lemma  8  Let  each  partition  Vh  consist  of  convex  polygons  {Ki],  i  —  1,  2, . . . ,  N, 
and  each  element  Ki  have  Nlc  corners  with  angles  j  =  1,2,...,  iV®  ( see  Fig¬ 
ure  3).  Let  wq  be  the  solution  to  (26)  for  a  given  q  G  H~a{Vh),  0  <  a  <  1. 

If  the  following  characteristic  equations  each  have  no  root  A  G  C  other  than 
A  =  —I  on  the  line  (a  —  2)1  (where  I  denotes  the  imaginary  variable): 

sinh2(AaA)  =  A2  sin2 (a;*),  j  —  1, . . . ,  Nlc,  %  —  1, . . . ,  N.  (28) 

Then,  every  wq  belongs  to  H4~a{Vh)-  Otherwise,  the  functions  wq  are  in 

H3~a(Vh). 


Proof:  This  lemma  is  a  result  of  Theorem  7. 2. 2. 3  and  Remark  7.2. 2. 4  in  the 
work  by  Grisvard  [19].  These  establish  the  regularity  of  solutions  to  the  bi¬ 
harmonic  equation  and  can  be  applied  to  establish  the  regularity  of  the  func- 


Fig.  3.  Corner  angles  of  the  polygons. 
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tions  wq,  as  these  satisfy  a  bi-harmonic  equation  on  each  element  Kf 

AAwq  =  $,  WKi  G  Vh,  (29) 

where  $  G  H~a(Ki )  and  homogeneous  Dirichlet  and  Neumann  conditions  hold 
on  d Ki.  To  prove  this  assertion,  we  start  by  substituting  test  functions  v  —  ip 
into  (26)  whose  restrictions  </y  to  Ji*  belong  to  V^Kf).  Thus,  we  get 

(A  (Pi,  A  Wq)  +  2  (Vipi,  Vwq)  +  (tfii,  Wq)  =  (fi,  q)ifi,  Wifi,  WKi  G  Vh., 

where  (•,  •)  denotes  the  duality  pairing  in  T>(Ki )  x  'D(Kj)' .  Application  of  the 
definition  of  the  distributional  derivative,  then  gives: 

(ifi,  AA  Wq  -  2  A  Wq  +  Wq)  =  (ifi,  (?)  ,  Wifi,  WKi  G  Vh. 

Comparison  with  (29)  reveals  that  <f>  =  q  +  2Awq  —  wq.  Since  q  G  H~a(Ki ) 
and  Awq,  wq  G  Lf^Kf),  <f>  must  belong  to  H~a(Ki).  The  proof  of  this  lemma  is 
completed  by  applying  Theorem  7.2. 2. 3  and  Remark  7. 2. 2. 4  in  [19]  to  (29)  and 
by  noting  that  the  functions  wq  satisfy  homogeneous  Dirichlet  and  Neumann 
boundary  conditions  on  dKt  and  that  every  polygon  Ki  is  convex  (he.  no  re¬ 
entrant  corners) .  ■ 


We  first  use  the  result  of  Lemma  8  to  derive  h  and  p  convergence  rates  in  the 
broken  space  Hl{Vh)- 

Theorem  4  (Convergence  in  Hl((Ph ) )  Let  u  G  H2((Ph)  be  the  unique  solution 
to  the  variational  problem  (8)  and  {uh  G  Vhp}  be  a  sequence  of  approxima¬ 
tions  (19)  using  families  of  regular  partitions  {Ph}  of  convex  polygons.  Then, 
the  approximation  error  eh  =  u  —  Uh  is  bounded  in  the  Hl(Ph)  norm  (25),  i.e. 


hi1  1  /  ~  ~ 

3C  >  0  :  \\eh\\H\vh)  A  C  IMI rrq/qp 

l>  V  Kievh 


P>  2,  n  >  2, 


where  r  =  min  {rA,  p  =  max  {pA,  h  =  max  {hA,  and  n  =  rnirhp  +  l,r). 

Ki evhl  5’  y  Ki£Phl  S  Ki&ThX  5’  P  v  ’ 

Proof:  We  return  to  (27)  and  set  a  =  1: 


\eh\\m(vh)  -  C 


h^v~ 4  /  ^  - 

VsW  L,  ll“lliri(jr,)  6y\P  , 

t  V  Kievh  ieH  Th) 


I  \\Wq\\2HH(Ki) 

Ki&Vh 


(30) 


We  now  call  upon  Lemma  8  to  determine  the  regularity  of  the  functions  wq, 
which  means  that  we  have  to  find  the  roots  A  G  C  to  the  following  character- 
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istic  equations  on  the  line  —I: 


sinh(Au;*)  =  ±A  sin(o;]),  j  —  1, . . . ,  Nlc,  i  —  1, . . . ,  N. 

Lemma  7.3. 2. 4  in  [19]  reveals  that  the  above  equations  have  only  one  root 
A  =  —  I  on  the  line  —  I.  From  Lemma  8  we  then  conclude  that  the  functions 
wq  are  in  H3{Vh)-  This  implies  that  s  >  3  and,  with  p  >  2,  we  obtain  v  >  3. 
Returning  to  (30),  this  gives: 


^  ^ 

eh\ \HHPh)  <  C—p  Y  IMlWd  suP 

P  V  Ki£Vh  4£H  1(ph> 


wq\\HHvh) 

IMI  H-HVh)’ 


(31) 


To  bound  the  term  involving  the  supremum,  we  return  to  (26)  and  rewrite 
this  as  follows: 


(v,  M*wq)  =  (v,  q),  Vu  G  H2(Vh ), 


where  B*  :  Hf){)(Vh)  — >  H~2(Vh )  is  a  linear  operator  associated  with  the 
bilinear  form  Since  we  know  that  the  bilinear  form  is  continuous  on 

H2(Vh)  x  H^0(Vh),  the  operator  B*(-)  is  continuous  on  H‘2}i){Vh)  r\H3(Vh),  he. 


\W*'Wq\\H-‘\Vh)  <  C  \\wq\\H2(Vh)  <  C  \\wq\\H3(ph). 

We  have  established  that  for  every  q  €  H^1  (Vh)  there  exists  a  unique  wq  G 
H3(Vh)nH20(Vh).  Hence,  !*(•)  is  bijective  from  H3(Vh)nH20(Vh)  to  H~\Vh) n 
H~2(Vh)-  Since  the  restriction  of  B*(-)  to  H3{Th)  is  both  continuous  and  bi¬ 
jective,  the  Banach  Theorem  states  that  the  inverse  of  the  restriction  of  B*(-) 
is  continuous.  Thus,  there  exists  C  >  0,  independent  of  q(-),  such  that: 

\\wq\\HHj>h)  =  II®"  ^lli?3^)  <  C  \\q\\H-*(vh),  ^ wq  ^  H3(Vh)  DH^Vh). 

By  backsubstituting  this  result  into  (31),  we  conclude  the  proof.  ■ 


Thus,  we  have  proved  optimal  h  and  p  convergence  of  the  error  in  i/1('P/l),  for 
p  >  2,  if  convex  polygonal  elements  are  used.  Optimal  h  and  p  convergence  in 
L2(fl)  is  also  possible  but  it  depends  on  the  shape  of  the  polygons  in  the  par¬ 
tition  Vh,  in  particular  the  value  of  the  corner  angles.  The  following  theorem 
is  a  consequence  of  Lemma  8  and  states  a  necessary  condition  on  the  corners 
of  the  polygons  in  order  to  obtain  optimal  convergence  rates  in  L2(h2). 

Theorem  5  (Convergence  in  L2( fi))  Let  u  G  H2(Vh )  be  the  unique  solution 
to  the  variational  problem  (8)  and  {uh  G  Vhp }  be  a  sequence  of  approxima¬ 
tions  (19)  using  families  of  regular  partitions  {Vh}  of  convex  polygons.  Let  Nlc 
denote  the  number  of  corners  of  an  element  Ki  and  denote  the  jth  corner 
angle  of  this  element.  If  the  following  characteristic  equations  have  no  roots  a 
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i»K  \  {0}, 

tanh(aa;*)  =  -  tan(2a;!)  a,  j  —  1,2, ,  Nlc,  i  =  1,  2, . . . ,  N.  (32) 
J  2  J 

Then,  the  approximation  error  eh  =  u  —  Uh  is  bounded  in  the  L2(f2)  norm  such 
that  3 C  >  0: 


II  II  L2(f2)  <  C 


h v 


-l 


P 

hP 


I  E  IMI^^y 

I<i£Vh 


e/i|U2(Q)  <  C  I  E  \\U\\ Hri{Ki)i 

p  V  Ki£Ph 


P  =  2,  Ti  >  2, 
P  >  3,  rt  >  2, 


where  r  =  min  {r^},  p  =  max  (p,; } ,  h  =  max{/q}7  and  fi  =  min(p  +  l,r). 

Kitzi'Ph 

Otherwise ,  the  error  is  bounded  by  the  following  suboptimal  bounds: 


h ^  i 

M U2(n)  <  c  ^n. 


E  IK2 


ATien 


P  >  2,  r *  >  2, 


Proof:  Since  ||-||L2(n)  <  ll'll-f/1^)’  the  error  in  L2(Tt)  converges  at  least  at  the 
same  rate  as  in  H1(Vh)-  To  obtain  a  lift  in  convergence  rates,  we  follow  the 
proof  of  Theorem  4,  recall  (27),  and  set  a  —  0,  i.e.: 


e/)||L2(n) 


C 


h^+ 


V- 4 


,r+s— 4 


P 


E  INI 

Ki&Vh 


2 

mi(Ki) 


sup 

q£L2(Q.) 


||  2 


lklU2(n) 


The  term  involving  the  supremum  can  be  bounded  in  the  same  manner  as 
previously  done  in  the  proof  of  Theorem  4,  but  we  again  have  to  determine 
the  regularity  of  the  functions  wq.  In  order  to  do  so,  we  call  upon  Lemma  8 
which  states  that  for  a  =  0  we  need  to  find  the  roots  A  G  C  on  the  line  —21 
of  the  following  characteristic  equations: 


sinh(Aa;j)  =  ±A  sin  (a;]),  j  =  1, . . . ,  Nlc,  i  —  1,. ...  ,N. 

By  expanding  A  =  a  —  21,  a  G  M,  we  get  sets  of  equations  that  govern  the  real 
and  imaginary  parts  of  the  above  equations, 


sinh(ao;*-)  cos (2a;*-)  =  ±asin(a;*), 
cosh(aa;*-)  sin(2a;*-)  =  ±2  sin  (a;]-). 


(33) 
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The  value  a  =  0  is  never  a  root  of  these  equations  for  0  <  a;*-  <  tt  (no  re-entrant 
corners),  as  substitution  of  a  =  0  into  (33) 2  gives: 

sin(2a/*)  =  ±2sin(o/*)  cos(cn))  =  ±1  a/*  =  0,7r. 

So  we  can  multiply  (33)2  by  a,  (33)1  by  2,  and  subtract  the  resulting  equa¬ 
tions,  which  yields: 

2  sinh(aa/*)  cos(2cn*)  =  ±a  cosh(aa/*)  sin (2a;*-),  a  0, 
which  we  can  rewrite  as  follows: 

tanh(ao/*)  =  -  tan(2a/*)  a,  a  0. 

Thus,  the  above  equations  are  equivalent  with  (28),  with  A  =  a  —  21.  If  there 
are  no  roots  a  other  than  a  =  0,  then  according  to  Lemma  8  the  functions  wq 
are  in  H4(Vh).  Thus,  for  p  —  2,  we  obtain  u  —  3,  and  for  p  >  3,  we  get  v  —  4, 
which  establishes  the  assertion.  ■ 

Corollary  6  Let  Vh  be  as  defined  in  Theorem  5.  In  addition,  assume  that  all 
corners  of  the  polygons  are  equal  to  n/2.  Then,  forp  >  3,  the  error  et  =  u  —  Uh 
converges  optimally  in  L2(Q),  i.e. 


where  r  =  min  in),  p  =  max  {pfit,  h  =  max  {hfit,  and  u  =  min(n  +  l,r). 

KitvhVlS,F  K^vhy- /yen1  J  P  u 


Proof:  This  is  an  immediate  consequence  of  Theorem  5,  as  substitution  of 
LUj  =  ti/ 2  into  (32)  yields: 


/  an  \ 

tanh  1  —  1  =0, 


which  has  only  one  root  a  =  0. 


From  Theorem  5  we  can  conclude  that  in  L2{fil)  we  always  obtain  suboptimal  h 
convergence  for  p  —  2.  For  p  >  3,  optimal  rates  can  be  obtained,  but  the  corner 
angles  have  to  satisfy  the  condition  (32).  This  condition  is  satisfied  for  strictly 
rectangular  meshes,  but  such  meshes  rarely  occur  in  practical  applications. 
Therefore,  only  in  few  occurrences  will  optimal  convergence  be  guaranteed  by 
Theorem  5. 
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4  Numerical  Verifications 

Jhl  One  Dimensional  Tests 


For  p  —  2,  3, 4,  5,  solutions  {ut}  of  (19)  are  computed  by  performing  successive 
uniform  h  refinements.  In  Figure  4,  the  convergence  results  are  shown  for  the 
norm  ||.||#2(-pft)  (5),  where  the  h  convergence  rates  are  computed  according  the 


following  rule: 


Ph 


log  (ej/e^1) 
log  2 


The  observed  convergence  rates  in  Figure  4  of  order  p  —  1  confirm  the  rates 
that  are  predicted  for  convergence  in  H2(Vh )  (see  Theorem  3).  In  Figures  5 
and  6,  the  results  are  illustrated  for  the  Hl{Vh)  and  L2(fl)  norm,  respectively. 
The  h  convergence  rates  in  Hl(Vh)  are  of  order  p  and  agree  with  the  prediction 
in  Theorem  4.  In  L2(fl),  the  convergence  rates  are  of  order  p  +  1,  for  p  >  3, 
and  of  order  p  for  p  =  2.  These  rates  also  confirm  the  rates  predicted  in 
Theorem  5  (for  one-dimensional  versions  of  (26),  the  dual  solutions  wq  are 
always  in  H4(Vh)  for  q(-)  that  belong  to  L2(fl)). 
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Convergence  Rate  H2(Ph)-Norms 
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Fig.  4.  Norm  |||-|||  of  the  approximation  error  (left)  and  uniform  h  convergence  rates 
(right)  versus  number  of  degrees  of  freedom  -  ID  results. 


H1(Ph) norm 


Convergence  Rate  H1(Ph)-Norms 
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Fig.  5.  Hl{Vh)  norm  of  the  approximation  error  (left)  and  uniform  h  convergence 
rates  (right)  versus  number  of  degrees  of  freedom  -  ID  results. 
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L2(0mega)  norm 


Convergence  Rate  L2(Omega)-Norms 


Fig.  6.  L2(Q)  norm  of  the  approximation  error  (left)  and  uniform  h  convergence 
rates  (right)  versus  number  of  degrees  of  freedom  -  ID  results. 
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4-2  Two  Dimensional  Tests 


Next,  we  consider  the  two-dimensional  test  case  on  the  unit  square  17  =  (0, 1)  x 

(0,1): 


u(x)  =  sin(27nr)  sinh(\/l  +  An2y)  (34) 

For  p  =  2,  3, 4, 5,  solutions  {uh}  of  (19)  are  computed  by  performing  succes¬ 
sive  uniform  h  refinements.  In  Figures  7  and  8,  the  convergence  results  are 
shown  for  the  norms  |||.|j|  and  |H|.F/2('Pfe),  respectively.  Note  that  the  norm  |j|'||j 
represents  a  broken  Laplacian  norm  on  Vh,  whereas  IHItf2^)  uses  the  com¬ 
plete  local  Sobolev  norm  in  H2,  i.e.  it  includes  the  local  L 2  norm  of  the  cross 
derivatives  d2u/dxdy  (see  also  (5)).  The  results  for  the  norm  IHIh2^)  confirm 
the  predicted  h  convergence  rates  of  Theorem  3.  The  convergence  rates  in  the 
norm  |||'|||  are  higher  in  the  pre-asymptotic  range  but  converge  to  the  same 
rates  as  those  observed  for  the  norm  IHIh2^)- 

Figure  9  shows  the  convergence  of  the  error  in  the  H1(Vh )  norm.  The  predicted 
rates,  stated  in  Theorem  4,  are  confirmed  by  exhibiting  convergence  rates  of 
order  p. 

Since  the  mesh  consists  of  rectangular  elements,  Theorem  5  and  Corollary  6 
assert  that  the  convergence  rates  in  L2(fl)  should  be  of  order  2,  for  p  =  2,  and 
p  +  1,  for  p  >  3.  The  convergence  results  shown  in  Figure  5  agree  with  this 
assertion. 

We  consider  an  additional  test  problem  on  the  quadrilateral  domain  depicted 
in  Figure  11.  The  domain  partitions  {Vh}  are  performed  as  illustrated  in  this 
figure.  Thus,  the  skewness  of  the  elements  is  determined  by  the  angle  6.  The 
same  Poisson  equation  as  in  the  previous  example  is  solved  on  the  quadrilateral 
and  the  Dirichlct  boundary  conditions  are  applied  such  that  the  solution  is  as 
given  in  (34). 

In  Figure  12,  the  uniform  h  convergence  rates  are  presented  for  6  =  n/6.  For 
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Convergence  Rate  Norm  |||e 


Fig.  7.  Norm  |||-|||  of  the  approximation  error  (left)  and  uniform  h  convergence  rates 
(right)  versus  number  of  degrees  of  freedom  -  2D  results. 
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Fig.  8.  H2(Vh)  norm  of  the  approximation  error  (left)  and  uniform  h  convergence 
rates  (right)  versus  number  of  degrees  of  freedom  -  2D  results. 
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HI  (Ph)  Norm  ||e 


Convergence  Rate  H1(Ph)  Norm 


dots 


dots 


Fig.  9.  norm  of  the  approximation  error  (left)  and  uniform  h  convergence 

rates  (right)  versus  number  of  degrees  of  freedom  -  2D  results. 


L2(Omega)  Norm  ||e||  Convergence  Rate  L2(Omega)  Norm 


Fig.  10.  L2(D)  norm  of  the  approximation  error  (left)  and  uniform  h  convergence 
rates  (right)  versus  number  of  degrees  of  freedom  -  2D  results. 
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convergence  in  H2(Vh)  and  //1( Vh),  both  figures  show  no  noticeable  difference 
with  the  results  given  in  Figures  7  through  9.  As  expected,  the  skewness  of  the 
mesh  does  not  affect  the  convergence  of  the  error  in  these  norms.  However, 
the  assertion  of  Theorem  5  for  the  convergence  in  the  L2(h2)  norm  of  skewed 
meshes,  is  only  confirmed  in  part.  The  corner  angles  in  the  mesh  are  either  7t/6 
or  57t/6.  For  the  latter  of  these,  the  characteristic  equations  (32)  have  nonzero 
roots.  Thus,  the  convergence  rates  in  L2(h2)  should  be  suboptimal  according 
Theorem  5.  For  the  even  order  approximations,  the  suboptimal  rate  of  order  p 
is  indeed  observed,  but  for  odd  order  approximation  the  optimal  rate  of  p  +  1 
is  still  obtained.  This  even-odd  behavior  in  L2(Q)  has  also  been  observed  for 
the  DGM  introduced  by  Oden,  Babuska,  and  Baumann  [4] ,  but  in  their  results 
the  suboptimal  L2(Q)  convergence  also  emerges  for  rectangular  meshes. 

5  Concluding  Remarks 

A  new  DGM  is  introduced  for  the  two-dimensional  reaction-diffusion  problem 
with  prescribed  Neumann  and/or  Dirichlet  boundary  conditions.  The  DGM 
formulation  employs  local  second  order  derivatives,  satisfies  a  local  conser¬ 
vation  property,  and  the  corresponding  bilinear  form  satisfies  coercivity  and 
continuity  conditions  on  the  broken  Sobolev  space  H2{Vh)-  Due  to  the  coer¬ 
civity  property,  the  formulation  is  numerically  stable  and  does  not  require  any 
additional  penalization. 

We  have  derived  a  priori  error  estimates  that  show  that  optimal  h  and  p 
convergence  is  obtained  in  H2(Vh)-  If  the  mesh  consists  of  convex  polygons, 
then  we  can  also  prove  optimal  convergence  in  Hl(Th)-  These  assertions  are 
confirmed  by  one-  and  two-dimensional  experiments. 
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Fig.  12.  Uniform  h  convergence  rates  versus  number  of  degrees  of  freedom  -  2D 
results  -  Mesh  distortion  of  30°. 
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In  L2(Q),  we  always  obtain  a  suboptimal  h  convergence  rate  of  order  p  for 
p  =  2.  For  p  >  3,  we  can  prove  optimal  h  and  p  convergence  if  the  corner 
angles  of  the  polygons  have  no  roots  other  than  zero  for  the  characteristic 
equation  (32).  If  we  employ  strictly  rectangular  meshes,  then  this  condition 
is  satisfied  and  optimal  convergence  for  p  >  3  is  obtained  and  confirmed  by 
numerical  results.  In  practical  applications,  meshes  consist  of  elements  with 
various  shapes  and  the  condition  (32)  is  most  likely  not  satisfied.  In  those 
cases,  we  can  only  prove  suboptimal  convergence  rates.  Remarkably,  for  odd 
order  approximations  on  skewed  meshes,  numerical  results  show  that  optimal 
convergence  is  obtained,  nevertheless.  This  ‘odd-even’  behavior  of  the  error  in 
L2(Vt)  is  also  observed  in  other  DGM’s  ( e.g .  see  [7]). 
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